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Motion of a continuous fluid can be decomposed into an "incompressible" rearrangement, which 
preserves the volume of each infinitesimal fluid element, and a gradient map that transfers fluid 
elements in a way unaffected by any pressure or elasticity (the polar decomposition of Y. Brenier). 
The Euler equation describes a system whose kinematics is dominated by incompressible rearrange- 
ment. The opposite limit, in which the incompressible component is negligible, corresponds to the 
Zel'dovich approximation, a model of motion of self-gravitating fluid in cosmology. 

We present a method of approximate reconstruction of the large-scale proper motions of matter 
in the Universe from the present-day mass density fleld. The method is based on recovering the 
corresponding gradient transfer map. We discuss its algorithmics, tests of the method against mock 
cosmological catalogues, and its application to observational data, which result in tight constraints 
on the mean mass density Sim and age of the Universe. 

PACS numbers: 95.35. +d, 95.75.Pq, 98.62.Py, 98.65.Dx 
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I. INTRODUCTION 

In the spectrum of possible models of fluid motion, the 
Euler equation of incompressible fluid constitutes an ex- 
treme. As was shown by Y. Brenier [ll[2], any Lagrangian 
motion of fluid admits a polar factorization into a compo- 
sition of an "incompressible" rearrangement, which pre- 
serves the volume of each inflnitesimal fluid element, and 
an "absolutely compressible" transfer, which displaces 
fluid elements to their flnal locations prescribed by the 
gradient of a suitable convex potential, while expanding 
or contracting them in a way unaffected by any pressure 
or elasticity. Decomposing a fuid motion into a sequence 
of small time steps and factoring out the compressible 
transfer from the inertial fluid motion at each step yields 
a difference scheme for the incompressible Euler equation 

In this article we show how the opposite approach, in 
which only the compressible transfer is retained, can be 
applied to solving the problem of reconstructing peculiar 
motions and velocities of dark matter elements in cosmol- 
ogy. We also discuss algorithmics of this method, which 
gives an explicit discrete approximation to polar decom- 
position and can also be applied to model incompressible 
fluid as suggested in [3]. 

Recall that on scales from several to a few dozen of 
h^^ Mpc the large-scale structure of the Universe is pri- 
marily determined by distribution of the dark matter. 
This distribution can be described by the mass density 
field and by the large-scale component of the peculiar 
velocity field PS^, controlled by dark matter itself via 



gravitational interaction. The dark matter distribution 
is traced by galaxies, whose positions and luminosities 
are presently summarized in extensive surveys [HOE]. 
On large scales luminosities of galaxies allow to deter- 
mine their masses, from which the mass density of the 
dark matter environment can be estimated using well- 
established techniques [7l[8]. 

It is appropriate to consider the reconstruction of the 
field of peculiar velocities as part of the more complex 
problem of reconstructing the full dynamical history of 
a particular patch of the Universe. Several approximate 
methods have been proposed to this end, of which we 
mention here two. The Numerical Action Method, based 
on looking for minimum or saddle-point solutions for a 
variational principle involving motion of discrete galax- 
ies, was introduced by P. J. E. Peebles in the late 1980s 
[9]; its modern state is addressed in the present vol- 
ume by A. Nusser |15j . In this paper we concentrate on 
the Monge- Ampere-Kantorovich method introduced in 
[TU] (hereafter the MAK method), specifically highlight- 
ing the structural relationship between the MAK method 
and the variational approach to the Euler equation of in- 
compressible fluid [TT]. 

In Section |IT] the mathematical setting of the MAK re- 
construction is derived by application of the Zel'dovich 
approximation to a suitable variational formulation of 
dark matter dynamics, which leads to the Monge- 
Kantorovich mass transfer problem and the Monge- 
Ampere equation. In Section [III] we discuss algorithmics 
of solving the disretized Monge-Kantorovich problem, 
which gives as a byproduct an algorithm of polar decom- 
position for maps between discrete flnite point sets. In 
Section IV we show that the MAK method performs very 
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well when tested against direct numerical simulations of 
the cosmological evolution and review the recent appli- 
cations of the MAK method to real observational data. 
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which yielded new tight constraints on the value of the 
mean mass density of the Universe. A detailed treatment 
of implementation and testing the MAK method against 
A'^-body simulations is presented in the companion paper 
by G. Lavaux in the present volume [20^. The paper is 
finished with a discussion and conclusions. 



II. DYNAMICS OF COLD DARK MATTER 
AND THE ZEL'DOVICH APPROXIMATION 

The most widely accepted explanation of the large- 
scale structure seen in galaxy surveys is that it results 
from small primordial fluctuations that grew under grav- 
itational self-interaction of coUisionless cold dark mat- 
ter particles in an expanding universe (see, e.g., [121 and 
references therein). The relevant equations of motion 
are the Euler-Poisson equations written here for a flat, 
matter-dominated Einstein-de Sitter universe (for a more 
general case see, e.g., [13]): 



drV + {v ■ S/xV) 
drP + Va; • (pv) 



Zt 

0, 
1 



(1) 
(2) 
(3) 



Here v{x,t) denotes the velocity, p{x,t) denotes the 
density (normalized so that the background density is 
unity) and 4>{x,t) is a gravitational potential. All quan- 
tities are expressed in comoving spatial coordinates x and 
linear growth factor r, which is used as the time vari- 
able; in particular, v is the Lagrangian r-time derivative 
of the comoving coordinate of a fluid element. A non- 
technical explanation of the meaning of these variables 
and a derivation of eqs. ([l]-|3| in the Newtonian approx- 
imation can be found, e.g., in [2]; see also [T5] in the 
present volume, where the growth factor is denoted by a. 

The right-hand sides of the Euler and Poisson equa- 
tions ([T]) and (|3| contain denominators proportional to 
T. Hence, it suffices for the problem not to be singular 
as T ^ that 

v{x,0) + W^cb{x,0) = 0, p{x,0) = l. (4) 

Note that the density contrast p — 1 vanishes initially, 
but the gravitational potential and the velocity, as de- 
flned here, stay finite thanks to our choice of the linear 
growth factor as the time variable. Eq. ^ provides ini- 
tial conditions at t = 0; at the present time r = tq 
the density is prescribed by a galaxy survey as explained 
above: 



pix,To) = pa{x). 



(5) 



In parallel with the Euler equation of incompressible 
fiuid, eq. ^ can be considered as the Euler-Lagrange 
equation for a suitable action O [TB] : 



Ia = l J^" dT J dx- t"(p|-(;|2 + a|V, 



(6) 



where a = | for the fiat Universe and minimization is 
performed under the constraints expressed by eqs. ([2||5|. 
Note that the term containing |Va;^p may be seen as a 
penalization for the nonuniformity of the mass distribu- 
tion, which corresponds to the lack of incompressibility 
of the fiuid; enhancing this penalization infinitely would 
suppress the "absolutely compressible" transfer of fiuid 
elements, thus recovering the incompressible Euler equa- 
tion. 

However, according to eq. Q the rotational component 
of the initial velocity field vanishes, which strongly sup- 
presses the "incompressible" mode of the fiuid motion at 
early times. Based on this observation, Ya. B. Zel'dovich 
proposed [17] an opposite approximation in which a 0. 
In this approximation eq. ([T]) assumes the form 



drV + {v ■ S/x)v = 0. 



(7) 



Much as the incompressible Euler system, the study 
of the Zel'dovich approximation benefits from the La- 
grangian approach. Let x{q,T) be the comoving coordi- 
nate at time r of a fiuid particle that was initially located 
at q: x{q,0) = q. Then 



p{x{q,T),T) 
v{x{q,T),T) 



{det{dx/dq)y 
drx{q,T), 



(8) 
(9) 



where the r derivative is taken as q is fixed. As observed 
by Zel'dovich, in these new variables the nonlinear eq. 
([t]) assumes a linear form 



dix = 0. 



(10) 



Moreover eq. ^ is satisfied automatically, and the action 
becomes 

Io=l J drjdq \drx{q, ^)\^ = ^ J dq \Mq) - 

(11) 

Here we denote xo{q) = x{q, tq) and use the fact that ac- 
tion minimizing trajectories of fiuid elements, determined 
by eq. (10 1, are given by 



x{q, T)^q+ {T/To){xo{q) - q) 



(12) 



Note that according to the first condition v{q,0) = 
{l/To){xo{q) — q) — V q^{q) and the Lagrangian map 
(12 1 remains curl-free for all t > 0: x{q,T) — q + 
T^gHq) = Vq$(q,T) with $(q,r) = \q\y2 + r^q); 
thus the "incompressible" rotational component of the 
fiuid motion is indeed fully suppressed. 

To find the motion of the fiuid in the Zel'dovich ap- 
proximation it is necessary to minimize the action (111 



under the constraint provided by the representation of 
density ([s]) and the boundary conditions (|4] |5]) : 



det{dxQ{q)/dq) = l/p{xo{q)). 



(13) 



In optimization theory this problem is called the Monge- 
Kantorovich problem. Equivalently, one can solve the 
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Monge- Ampere equation that follows from (13 1 for the 
function <i>o('7) = ^{QiTq) such that Xo{q) = Vq<f>o(9): 

det{dHo{q) /^q^^qJ) = l/po(Vq<i>o(g)). (14) 

At large scales the Lagrangian map Xo{q) in cosmology is 
free from multistreaming (the presence of several streams 
of dark matter at the same spatial location) . Under this 
assumption the potential <&o('7) is necessarily convex [37], 
and the Legendre transform 



^'o(a;) = max(q • x - $0(9)), 



(15) 



where the maximum is attained at q such that x 
Vq<l>o(<7), gives (14) a simpler form 



det{d'^'l!o{x)/dx,dx.j) ^ po{x). 



(16) 



The Monge- Ampere-Kantorovich (MAK) method, intro- 
duced in [TU] , consists in solving either of these two prob- 
lems for a;o(q) and using eqs. (12] 9| to approximately 
recover the present field v{x,Tajoi peculiar velocities. 



III. ALGORITHMICS OF SOLVING 
THE MONGE-KANTOROVICH PROBLEM AND 
THE DISCRETE POLAR DECOMPOSITION 

To solve the Monge-Kantorovich problem numerically 
we discretize the initial and final distributions of mass 
into collections of Dirac point masses: all inital point 
masses (/i, qj) are assumed to lie on a regular grid and 
be equal, whereas the masses {mj,Xj) discretizing the 
present distribution po{x) typically come from a galaxy 
survey. The discretized action functional (111 assumes 
the form 



(17) 



where 7,^- > denotes the amount of mass transferred 
from qi to Xj and mass conservation implies for all i,j 



J2 



p.. 



(18) 



In practice we choose all rrij to be integer multiples of the 
elementary mass /i, which guarantees that all 7^^ assume 
only values or /i. 

Observe that the unknowns 7^^ enter into the problem 
of minimizing (17 1 under constraints (18 1 linearly. Prob- 



lems of this form are called linear programs and can be 
efficiently solved by various optimization methods, see, 
6-g-, |18l . Often the original, primal formulation of a lin- 
ear program is treated simultaneously with a (Lagrange) 
dual formulation, which is another linear program; such 
algorithms are called primal-dual algorithms. For the Hn- 
ear program at hand the dual formulation is to maximize 



(19) 



under the constraints 
1 



\xj — qi\^ + (pi + ^pj > for all 



(20) 



Here <f>i, 4'j are Lagrange multipliers for constraints (18 1 



the duality comes from the following representation of 
the (coinciding) optimal values of the two problems: 



id 



mm max 

7ij>o 0i.V'3 \ 2 



where taking max or min leads to the primal and dual 
problem respectively. The coincidence of the optimal val- 
ues implies that 7^-, 4>i, ipj solve the respective problems 
if and only if 



(21) 



In view of the nonnegativity conditions this means that 
for each pair {i,j) either jij = or constraint (20 1 is 
satisfied with equality (and then jij = fi). For all other 
values of 7^^ > and (/)i,^^that satisfy constraints (18 
20 1, the left-hand side of (21 1 is strictly positive and thus 



the value of ( 17) is strictly greater than that of ( 19 1; such 
7ij, Tji, ipj cannot be solutions to the respective optimiza- 
tion problems. 

A typical primal-dual algorihtm starts with a set of 
values of (pi , ipj for which all inequalities ( 20 1 hold and 



proceeds in a series of steps. At each step a constraint of 
the form (20) is found that is, in a certain sense, the eas- 



iest to be turned into equality, values of the correspond- 
ing (pi and are updated accordingly, and the 7^ is set 
to p. An algorithm stops when the number of equalities 



in (20) equals the number of masses (/Lt, qi), so that (21 1 
is satisfied. 

We found the auction algorithm of D. Bertsekas [TS] 
(see also [2^) to be a particularly efficient primal-dual 
algorithm for the huge data sets arising from the cosmo- 
logical application. The search for the constraint (20 1 



that is to be satisfied with equality at each step may be 
performed very efficiently using a specially developed ge- 
ometrical search routine, which is based on suitably mod- 
ified routines of the ANN library [2T] . Further details of 
this numerical implementation of the MAK method are 
given in [14] and in our forthcoming publication with 
M. Henon. A new implementation in a parallel comput- 
ing environment is reported in the companion paper of 
G. Lavaux [20 . 

We finally show why the solution of the Monge- 
Kantorovich problem in the discrete case gives a discrete 
analogue of polar decomposition. Let a discrete "map" 
7* between two sets of points (/i, Qi) and {mj,Xj) be 
given such that rrij — J2i llj ^^id 7,*- take only values 
and /i. Solving the corresponding Monge-Kantorovich 
problem will give a discrete analogue 7ij of the gradient 
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TABLE I: MAK reconstruction: percentage of successfully 
reconstructed initial positions at different scales (measured in 
units of mesh size in a box of 128^ grid points). 



logio (1 + surface density ot pcdnts [ erbitrery ]) 



Scale 


0" 


< 1 


< 2 


< 3 


< 4 


< 5 


% 


18% 


41% 


54% 


66% 


74% 


81% 



"Exact reconstruction. 



transfer. To see this observe that for <I>i 



^3 



\\Xj\'^ 



tpj equahty in ( 20 1 means that 



l\q^ 



max(<7fc • Xj 

k 



$fc) 



(22) 



(cf. (lisp), i.e., the map sending qj to X j IS cl discrete ana- 
logueoi the gradient map x{q) — Vq<f>(q) of the previous 
section. The corresponding discrete analogue of "incom- 
pressible" rearrangement is a permutation of (/i, q^) that 
for any j sends the set of points i such that 7*^ > to the 
set of points i' such that ji'j > 0; if furthermore nij = fj, 
for all j, both sets are singletons and the permutation is 
recovered uniquely. 

In the MAK method we are interested in the "gradi- 
ent" part of the Lagrangian map sending elements of dark 
matter to their present positions, whereas in the differ- 
ence scheme for the Euler equation proposed by Y. Bre- 
nier in [7 it is the permutation part that is retained. 

Finally note that an alternative approaches to solve the 
Monge-Kantorovich problem, based on direct numerical 
resolution of eqs. (l7l[2|), was proposed in 



IV. 



TESTING AND APPLICATION OF THE 
MAK METHOD TO COSMOLOGICAL 
RECONSTRUCTION 



The validity of the MAK method depends on the qual- 
ity of the Zel'dovich approximation, which is hard to es- 
tablish rigorously. To be able to apply the method to 
real-world data we have instead to rely on extensive nu- 
merical tests. 

We report here a test of the MAK method against an 
iV-body simulation with over 2- 10^ particles The A^- 
body simulation had the following characteristics: 128'^ 
particles were assembled in a cubic box of 200 Mpc, 
giving the mean inter-particle separation of 1.5 Mpc; 
the initial conditions for the velocity field v{q,0) were 
taken Gaussian; the density parameter was chosen to be 
il,n = 0.3, the Hubble parameter to be /i = 0.65, the 
normalization of the initial power-spectrum erg = 0.99, 
and the mass of a single particle in this simulation was 
M = 3.2 X IO^^H'^Mq where Mq is the solar mass. 

The scatter plot in Fig. [T] demonstrates the perfor- 
mance of the MAK method in finding Lagrangian posi- 
tions of the particles. At the scales that were probed, po- 
sitions of about 20% of particles are reconstructed exactly 
(for detailed data see Table This low rate is due to 
large multi-streaming at such small scales; at larger scales 
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FIG. 1: Scatter plot of the MAK-reconstructed initial coor- 
dinates of particles versus their true initial coordinates for 
a sample of an A-body simulation with 128"^ particles in a 
cubic box of size 200 Mpc; ideal reconstruction would 
correspond to the diagonal. A "quasi-periodic (QP) projec- 
tion" coordinate q = {qi + g2\/2 -I- q:i\/3)/{l + y/2 + y/S) is 
used with < qi < 1, where 1 corresponds to the (rescaled) 
box size. The QP projection maps a regular grid in the unit 
cube into the unit segment such that q images of no two grid 
points coincide. Shown is the decimal logarithm of the local 
density of points plus 1; the resolution in QP coordinates is 
1/256, the Lagrangian mesh spacing is 1/128. 



where mean inter-particle separation is about 6h~^ Mpc 
(up to 3 meshes), the MAK reconstruction gives the La- 
grangian positions of two thirds of the particles exactly. 

At the moment of its implementation in 2006 this re- 
construction of 128'^ particles was unprecedented and 
broke a computational barrier for cosmological recon- 
struction schemes; more detailed tests of the MAK re- 
construction of the same scale are reported by G. Lavaux 
pO in the present volume. 

We now turn to applications of the MAK reconstruc- 
tion to real observational data. The MAK reconstruction 
of the peculiar velocities depends on the mean density of 
the Universe (the parameter fim), the age of the Universe 
To (to ^ I/Hq where Hq = 100 xhis the Hubble parame- 
ter), and the mass-luminosity relation. Assuming certain 
values of these parameters, one can estimate the peculiar 
velocities of the galaxies as ratios of their displacements 
to tq. Optimizing the matching between these velocities 
and the observed velocities of a few test galaxies, one can 
then constrain the parameters of the reconstruction. 
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Age (Gyrs) Method 

FIG. 2: Constraints on the mean mass density and the age of 
the Universe obtained by applying reconstruction techniques 
to real observational data. Left pane: solid contours mark o 
and 2(T confidence levels for the MAK reconstruction (shaded) 
and the least-action (LA) reconstruction (unfilled). Shaded 
also is the confluence of the constraints on density and age pa- 
rameters from WMAP [24| of Qmh^ = 0.134 and from SDSS 
[B] of Q.mh — 0.21. The 2a concordance region of the four 
methods is filled. Right pane: density parameter estimates 
for Ho = 80. For details see [25]. 

This procedure is illustrated in Fig. [2] taken from [25] . 
The catalog of galaxies {rrij , Xj ) that is used here is a 40% 
augmentation of the Nearby Galaxies Catalog, now in- 
cluding 3300 galaxies within 3000 km s~^ [1]. This depth 
is more than twice the distance of the dominant com- 
ponent, the Virgo cluster, and the completion to this 
depth in the current catalog compares favourably with 
other all-sky surveys. The second observational com- 
ponent is an extended catalog of galaxy distances (or 
radial component of peculiar velocities). In this cata- 
log, there are over 1400 galaxies with distance measures 
within the 3000 km s~^ volume; over 400 of these are de- 
rived by high quality observational techniques that give 



accurate estimates of the radial components of peculiar 
velocities. The important feature of Fig. |2] is that the 
MAK contours are transversal to contours provided by 
other methods, which largely reduces the degeneracy of 
constraints in the parameter space. 



V. CONCLUSION 

According to the polar decomposition theorem of 
Y. Brenier [Hd], kinematics of continuous fluid motion 
can be decomposed into "incompressible" rearrangement 
and "infinitely compressible" gradient transfer. The Eu- 
ler equation of describes a system whose kinematics is 
dominated by incompressible motion. In this paper we 
show that the opposite limit, in which the incompressible 
component is negligible, corresponds to the Zel'dovich 
approximation, a physically meaningful model of motion 
of self-gravitating fluid arising in cosmology. 

This result enables us to approximately reconstruct pe- 
culiar motion of matter elements in the Universe from 
information on their present-day distribution, without 
any knowledge of the velocity field; indeed, the latter it- 
self can be recovered from the reconstructed Lagrangian 
map. The viability of this method is established by test- 
ing it against a large-scale direct A^-body simulation of 
cosmological evolution; when applied to real obsrvational 
data, the method allows to get very tight constraints on 
the valuse of the mean mass density and gthe age of the 
Universe. 

Another our contribution is an efficient numerical 
method decomposing a given displacement field into "in- 
compressible" and "infinitely compressible" parts. This 
method is not limited to cosmological reconstruction but 
can also be used for modelling the dynamics of incom- 
pressible fluid as suggested by Y. Brenier in [3^ . 
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